Implementation of the generalised free volume theory
Fleer, G. J., & Tuinier, R. (2008). Analytical phase diagrams for colloids and non-adsorbing polymer. Advances in Colloid and Interface Science, 143(1-2), 1–47. doi:10.1016/j.cis.2008.07.001
In [1]:
from matplotlib.pyplot import *
%matplotlib inline
from colloids import phase
We use the colloids.phase module, where the main object represents an equation of state of the colloids. Here we will use the standard Carnahan and Starling equation of state for the hard sphere fluid and the Hall equation of state for the hard sphere crystal. These objects are almost free to construct so we will call the constructor explicitely every time.
To take the polymers into account, one need to know the size ratio $q_R = 2R_g/\sigma_C$ where $\sigma_C$ is the diameter of the colloids and $R_g$ is the radius of gyration of the polymers.
Calculations are done in a parameter space $(f, \Pi v)$ convinient for the theory and then converted to the experimental parameter space $(\phi, C_p)$.
First, we compute the fluid-solid coexistence for $q_R=0.1$ in the $(f, \Pi v)$ space. For each $\Pi v$ between 0 and maxpiv we finding the $f_F$ and $f_S$ that satisfy the equality of the chemical potentials and of the pressures.
In [2]:
qR = 0.10
q = phase.qR2q(qR)
FS = phase.all_coexistence(q, phase.CarnahanStarling(), phase.Hall(), maxpiv=1000)
#tie lines
for piv, fF, fS in FS[::10]:
plot([fF, fS], [piv]*2, color=[0.9]*3)
#coexistence lines
for f in FS.T[1:]:
plot(f, FS[:,0], color='k')
xlabel(r'$f$')
ylabel(r'$\Pi v$')
Out[2]:
We can then convert $f$ to a more readily understandable volume fraction $\phi$.
In [3]:
#tie lines
for piv, fF, fS in FS[::10]:
plot([phase.f2vf(fF), phase.f2vf(fS)], [piv]*2, color=[0.9]*3)
#coexistence lines
for f in FS.T[1:]:
plot(phase.f2vf(f), FS[:,0], color='k')
xlabel(r'$\phi$')
ylabel(r'$\Pi v$')
Out[3]:
The conversion from $\Pi v$ to polymer concentration $C_P$ is not that straightforward since we also need the value of $f$. And we only obtain the $C_P/C_P^*$ where $C_P^*$ is the polymer concentration at overlap.
Note that due to the dependence on $f$, the tie lines now have a slope.
In [4]:
#tie lines
for piv, fF, fS in FS[::10]:
plot(
[phase.f2vf(fF), phase.f2vf(fS)],
[phase.piv2y(piv, qR) * phase.alpha(fF, q), phase.piv2y(piv, qR) * phase.alpha(fS, q)],
color=[0.9]*3
)
#coexistence lines
for f in FS.T[1:]:
plot(phase.f2vf(f), phase.piv2y(FS[:,0], qR) * phase.alpha(f, q), color='k')
xlabel(r'$\phi$')
ylabel(r'$C_P/C_P^*$')
Out[4]:
Of course, this can be done for different size ratios, but we have to be careful about the value of maxpiv. If too small, the algorithm may not converge.
In [5]:
figure(figsize=(16,6))
for i,qR in enumerate([0.04,0.06,0.08, 0.1,0.2]):
q = phase.qR2q(qR)
FS = phase.all_coexistence(q, phase.CarnahanStarling(), phase.Hall(), maxpiv=1./qR**3)
#coexistence lines
for f in FS.T[1:]:
subplot(1,2,1).plot(f, FS[:,0], color=cm.autumn(i/5.))
subplot(1,2,2).plot(phase.f2vf(f), phase.piv2y(FS[:,0], qR) * phase.alpha(f, q), color=cm.autumn(i/5.))
subplot(1,2,1)
ylim(0,1000)
xlabel(r'$f$')
ylabel(r'$\Pi v$')
subplot(1,2,2)
ylim(0,0.5)
xlabel(r'$\phi$')
ylabel(r'$C_P/C_P^*$')
Out[5]:
Given sufficient attraction between the particles, here osmotic pressure exerted by the polymers, a low density fluid (gas) can cohexist with a high density fluid (liquid). This time we explore intrinsic properties of the fluid equation of state, thus we will use only Carnahan and Starling.
We obtain the location of the critical point in the $(f, \Pi v)$ space and convert it as before.
In [6]:
qR = 0.15
q = phase.qR2q(qR)
fc, pivc = phase.CarnahanStarling().critical_point(q)
print fc, pivc
print phase.f2vf(fc), phase.piv2y(pivc, qR) * phase.alpha(fc, q)
For each $\Pi v$ between the critical pressure and maxpiv we find the $f_G$ and $f_L$ that satisfy the equality of the chemical potentials and of the pressures. We obtain also the location of the spinodal lines as respective location of minima and extrema of the free energy.
In [7]:
figure(figsize=(16,6))
GL = phase.CarnahanStarling().all_GL(q, 200)
#tie lines
for piv, fG, fL in GL[::10,:3]:
subplot(1,2,1).plot([fG, fL], [piv]*2, color=[0.9]*3)
subplot(1,2,2).plot(
[phase.f2vf(fG), phase.f2vf(fL)],
[phase.piv2y(piv, qR) * phase.alpha(fG, q), phase.piv2y(piv, qR) * phase.alpha(fL, q)],
color=[0.9]*3
)
#binodal lines
for f in GL.T[1:3]:
subplot(1,2,1).plot(f, GL[:,0], color='k')
subplot(1,2,2).plot(phase.f2vf(f), phase.piv2y(GL[:,0], qR) * phase.alpha(f, q), color='k')
#spinodal lines
for f in GL.T[3:]:
subplot(1,2,1).plot(f, GL[:,0], 'k--')
subplot(1,2,2).plot(phase.f2vf(f), phase.piv2y(GL[:,0], qR) * phase.alpha(f, q), 'k--')
subplot(1,2,1)
ylim(0,200)
axvline(phase.vf2f(0.64))
xlabel(r'$f$')
ylabel(r'$\Pi v$')
subplot(1,2,2)
ylim(0,0.6)
axvline(0.64)
xlabel(r'$\phi$')
ylabel(r'$C_P/C_P^*$')
Out[7]:
Do not look too much at large volume fractions: the liquid branch of the binodal and the spinodal are going to volume fractions over close packing (vertical blue line)! This is obviously unphysical and show the limit of Carnahan and Starling equation of state that has no divergence at random close packing.
Liu equation of state is much better behaved. Beware, long evaluation
In [8]:
figure(figsize=(16,6))
for c, eos in zip([np.zeros(3), np.array([1,0,0])], [phase.CarnahanStarling(), phase.Liu()]):
GL = eos.all_GL(q, 200)
#tie lines
for piv, fG, fL in GL[::10,:3]:
subplot(1,2,1).plot([fG, fL], [piv]*2, color=np.maximum(c,0.9))
subplot(1,2,2).plot(
[phase.f2vf(fG), phase.f2vf(fL)],
[phase.piv2y(piv, qR) * phase.alpha(fG, q), phase.piv2y(piv, qR) * phase.alpha(fL, q)],
color=np.maximum(c,0.9)
)
#binodal lines
for f in GL.T[1:3]:
subplot(1,2,1).plot(f, GL[:,0], color=c)
subplot(1,2,2).plot(phase.f2vf(f), phase.piv2y(GL[:,0], qR) * phase.alpha(f, q), color=c)
#spinodal lines
for f in GL.T[3:]:
subplot(1,2,1).plot(f, GL[:,0], color=c, ls='--')
subplot(1,2,2).plot(phase.f2vf(f), phase.piv2y(GL[:,0], qR) * phase.alpha(f, q), color=c, ls='--')
subplot(1,2,1)
ylim(0,200)
xlabel(r'$f$')
ylabel(r'$\Pi v$')
axvline(phase.vf2f(0.64))
subplot(1,2,2)
ylim(0,0.6)
axvline(0.64)
xlabel(r'$\phi$')
ylabel(r'$C_P/C_P^*$')
Out[8]:
In [16]:
for i,qR in enumerate([0.071, 0.087, 0.1, 0.15]):
q = phase.qR2q(qR)
GL = phase.CarnahanStarling().all_GL(q, 1./qR**3)
for f in GL.T[3:4]:
plot(phase.f2vf(f), phase.piv2y(GL[:,0], qR) * phase.alpha(f, q), color=cm.autumn(i/4.))
ylim(0,0.6)
xlabel(r'$\phi$')
ylabel(r'$C_P/C_P^*$')
Out[16]:
Here we have two cases, depending on $q_R$. For small $q_R$, the critical point is over the fluid line, thus the gas-liquid coexistence is metastable to crystallisation.
In [65]:
qR = 0.15
q = phase.qR2q(qR)
figure(figsize=(16,6))
GL = phase.CarnahanStarling().all_GL(q, 200)
#binodal lines
for f in GL.T[1:3]:
subplot(1,2,1).plot(f, GL[:,0], color='k')
subplot(1,2,2).plot(phase.f2vf(f), phase.piv2y(GL[:,0], qR) * phase.alpha(f, q), color='k')
#spinodal lines
for f in GL.T[3:]:
subplot(1,2,1).plot(f, GL[:,0], 'k--')
subplot(1,2,2).plot(phase.f2vf(f), phase.piv2y(GL[:,0], qR) * phase.alpha(f, q), 'k--')
FS = phase.all_coexistence(q, phase.CarnahanStarling(), phase.Hall(), maxpiv=200)
#coexistence lines
for f in FS.T[1:]:
subplot(1,2,1).plot(f, FS[:,0], color='k')
subplot(1,2,2).plot(phase.f2vf(f), phase.piv2y(FS[:,0], qR) * phase.alpha(f, q), color='k')
subplot(1,2,1)
ylim(0,200)
axvline(phase.vf2f(0.64))
xlabel(r'$f$')
ylabel(r'$\Pi v$')
subplot(1,2,2)
ylim(0,0.6)
axvline(0.64)
xlabel(r'$\phi$')
ylabel(r'$C_P/C_P^*$')
Out[65]:
Or, with Liu equartion of state
In [67]:
qR = 0.15
q = phase.qR2q(qR)
figure(figsize=(16,6))
GL = phase.Liu().all_GL(q, 200)
#binodal lines
for f in GL.T[1:3]:
subplot(1,2,1).plot(f, GL[:,0], color='k')
subplot(1,2,2).plot(phase.f2vf(f), phase.piv2y(GL[:,0], qR) * phase.alpha(f, q), color='k')
#spinodal lines
for f in GL.T[3:]:
subplot(1,2,1).plot(f, GL[:,0], 'k--')
subplot(1,2,2).plot(phase.f2vf(f), phase.piv2y(GL[:,0], qR) * phase.alpha(f, q), 'k--')
FS = phase.all_coexistence(q, phase.Liu(), phase.Hall(), maxpiv=200)
#coexistence lines
for f in FS.T[1:]:
subplot(1,2,1).plot(f, FS[:,0], color='k')
subplot(1,2,2).plot(phase.f2vf(f), phase.piv2y(FS[:,0], qR) * phase.alpha(f, q), color='k')
subplot(1,2,1)
ylim(0,200)
axvline(phase.vf2f(0.64))
xlabel(r'$f$')
ylabel(r'$\Pi v$')
subplot(1,2,2)
ylim(0,0.6)
axvline(0.64)
xlabel(r'$\phi$')
ylabel(r'$C_P/C_P^*$')
Out[67]:
For larger $q_R$ the critical point is below the fluid line, so we can have triple coexistence between gas, liquid and fluid
In [96]:
qR = 0.3
q = phase.qR2q(qR)
figure(figsize=(16,6))
GL = phase.CarnahanStarling().all_GL(q, 50)
#binodal lines
for f in GL.T[1:3]:
subplot(1,2,1).plot(f, GL[:,0], color='k')
subplot(1,2,2).plot(phase.f2vf(f), phase.piv2y(GL[:,0], qR) * phase.alpha(f, q), color='k')
#spinodal lines
for f in GL.T[3:]:
subplot(1,2,1).plot(f, GL[:,0], 'k--')
subplot(1,2,2).plot(phase.f2vf(f), phase.piv2y(GL[:,0], qR) * phase.alpha(f, q), 'k--')
FS = phase.all_coexistence(q, phase.CarnahanStarling(), phase.Hall(), maxpiv=50)
#coexistence lines
for f in FS.T[1:]:
subplot(1,2,1).plot(f, FS[:,0], color='k')
subplot(1,2,2).plot(phase.f2vf(f), phase.piv2y(FS[:,0], qR) * phase.alpha(f, q), color='k')
subplot(1,2,1)
ylim(0,50)
axvline(phase.vf2f(0.64))
axhline(phase.CarnahanStarling().critical_point(q)[1]*1.1)
xlabel(r'$f$')
ylabel(r'$\Pi v$')
subplot(1,2,2)
ylim(0,0.6)
axvline(0.64)
xlabel(r'$\phi$')
ylabel(r'$C_P/C_P^*$')
Out[96]:
TO DO: fix instability problem of the GL near critical point for $q_R>0.2$. Implement triple point
In [86]:
reload(phase)
Out[86]:
In [18]:
qR = 0.088
q = phase.qR2q(qR)
print q
In [19]:
print phase.CarnahanStarling().critical_point(q)
In [20]:
#Compute and plot spinodal
GL = phase.CarnahanStarling().all_GL(q, 3000)
plot(phase.f2vf(GL[:,3]), phase.piv2y(GL[:,0], qR)*phase.alpha(GL[:,3],q), 'k-')
#samples
for name, c in zip(['gel', 'fluid', 'clusters'], ['rs', 'bo', 'yd']):
phi, cp = np.loadtxt(
'/home/mathieu/Documents/Recherches/tex/gel_cond/trunk/phase_diag_{}.csv'.format(name),
usecols=[1,2], skiprows=1,
unpack=True)
plot(phi/100., cp, c)
#sgp = np.loadtxt('/home/mathieu/Documents/Recherches/tex/gel_cond/trunk/gasliquid_sg.phd')
#plot(sgp[:,0], sgp[:,1], 'y--')
xlim(0,0.4)
ylim(0,2)
Out[20]:
In [50]:
(0.12/0.9)**(1/0.9)
Out[50]:
In [ ]: